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We discuss the relation between the fluctuation-dissipation relation derived by Chatelain and 
Ricci-Tersenghi [C. Chatelain, J.Phys. A 36, 10739 (2003); F. Ricci-Tersenghi, Phys.Rev.E 68, 
065104(R) (2003)] and that by Lippiello-Corberi-Zannetti [E. Lippiello, F. Corberi and M. Zannetti 
Phys. Rev. E 71, 036104 (2005)]. In order to do that, we re-derive the fluctuation-dissipation 
' relation for systems of discrete variables evolving in discrete time via a stochastic non-equilibrium 

^ , Markov process. The calculation is carried out in a general formalism comprising the Chatelain, 

■ Ricci-Tersenghi result and that by Lippiello-Corberi-Zannetti as special cases. The applicability, 

generality, and experimental feasibility of the two approaches is thoroughly discussed. Extending 
the analytical calculation to the variance of the response function we show the vantage of fleld-free 
. numerical methods with respect to the standard method where the perturbation is applied. We also 

show that the signal to noise ratio is better (by a factor \/2) in the algorithm of Lippiello-Corberi- 
Zannetti with respect to that of Chatelain- Ricci Tersenghi. 
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I. INTRODUCTION 



Recently, there has been much interest in the extension in the out of equilibrium regime of the fluctuation-dissipation 
, , theorem (FDT) , through more general fluctuation-dissipation relation (FDR) , which have led to the concept of effective 
^ ' temperature [l| and to the connection between non-equilibrium and equilibrium properties 15] . Fluctuating two-time 
, quantities have also been actively investigated, particularly in relation to the detection and quantification of dynamical 
' heterogeneities, mostly in disordered systems [j]- 

The search of FDR between response functions and properties of the unperturbed system, has led to a number of 
proposals Among these, the two by Chatelain [6] and Ricci-Tersenghi J] (CRT) and the one by Lippiello, 

Corberi, Zannetti [1, ^ij, [l^] (LCZ) have succeeded in making the connection between the dynamical susceptibility 
and unperturbed correlators between observable quantities. In addition to the intrinsic theoretical interest, these 
results opened the way to the development of perturbation-free numerical algorithms, allowing for highly efficient and 
precise measurements of the response function via correlators, without need of switching on any perturbation. 

However, the paths followed by CRT on one side and by LCZ on the other, are quite different as well as the 
final results. The two approaches lead to expressions of the susceptibility in terms of radically different unperturbed 
; ^ correlation functions, making the mapping between them cumbersome. This poses the question of understanding the 
inner relationship between the two results, of their degree of generality, of which is performing better in numerical 
implementations and of the possible experimental implications. In this paper we study these issues and answer these 
questions. In order to carry out this program, we derive the FDR for systems evolving in discrete time via stochastic 
Markov processes, defined by transition probabilities obeying detailed balance. We develop a unified formalism 
containing different approaches as special cases and explain the difference between those of CRT and LCZ: while 
in the LCZ case the response function is related to correlation functions computed over the whole non-equilibrium 
ensemble, in the CRT approach, instead, averages are taken over a restricted set of trajectories. 

The derivation is fully general for what concerns the nature of the discrete variables (e.g. Ising, Potts, Clock 
etc.) and of the transition probabilities. However, the constraint of a restricted set of trajectories in the CRT 
approach requires the microscopic knowledge of the sequence of (attempted) updates, which is manageable only in 
numerical simulations. On the other hand, in the LCZ approach a standard unrestricted ensemble average is involved, 
and the response function is written in terms of standard correlation functions between observable quantities. This 
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allows analytical treatments by means of the usual methods of statistical mechanics and, in principle, experimental 
applications. On the other hand, other approaches, such as those in |^, [lo| do not express the response function in 
terms of observables. 

After clarifying the relations between the FDR in the CRT and LCZ approaches, we turn to compare the efficiencies 
of the numerical algorithms based on them, together with that of the standard method (SM) , requiring the application 
of an external perturbation h. An important advantage of the perturbation-free methods is that the limit /i — > is 
built in, while, with the SM, checking for linearity is often numerically demanding. Besides this, field-free methods 
are also characterized by a better signal to noise ratio. This is a relevant fact, since the numerical computation of 
the response function is extremely noisy. In order to quantify such a noise, we compute exactly the variance of the 
response function for each of the three algorithms. In the SM it diverges as preventing small values of h to be 

used and making linearity often insecure. Since, obviously, such a drawback is absent with the field-free methods, 
there is an enormous advantage in their implementation. Nevertheless, also with perturbation-free algorithms the 
noise can be significant, especially for large time differences. The comparison between the variances of the CRT and 
the LCZ methods shows that the LCZ approach yields a better signal to noise ratio by a factor \/2. Basically, this 
difference is a consequence of the restriction in the set of trajectories required by the CRT method. 

In addition to the relevance for numerical applications, the results for the variances give a contribution to the 
understanding of fluctuations of two-time quantities, shedding some light in the field of nonlinear susceptibilities p^ . 

After investigating and clarifying the relation between the different algorithms and their performances, we present 
the results of numerical simulations in order to discuss the generality of the method and to illustrate the efficiency of 
the field-free algorithms with particular examples. We compute numerically the response function for models of Ising 
spins (the ferromagnetic Ising model and the Edwards- Anderson (EA) spin glass in d = 3) with the three methods. 
We show that both the CRT and the LCZ algorithms produce, with great accuracy, the same response which can be 
obtained with the SM. Computing the variances of the three methods, we obtain the results outlined above. Finally, 
we compute the response function in the Fredrickson- Andersen (FA) model, both applying the perturbation and 
with the field-free method of LCZ, finding again perfect agreement. This demonstrates the applicability of the LCZ 
algorithm also in this case, and that the criticism raised in Ref. fTJ| does not hold. 

This Article is organised as follows: In Sec.|n]we present the derivation of the FDR. We discuss the results obtained, 
their generality and the measurability of the correlators involved. In Sec. IIIII we compute and compare the variances 
of the three algorithms. Sec. llVl is devoted to explicit numerical implementations; We consider the 3c? ferromagnetic 
Ising and EA models quenched to the critical temperature and below it. Sec. II V B] contains the application to the FA 
model. The conclusions are drawn in Sec. |Vl where some perspectives are discussed. 

II. ANALYTICAL DERIVATION OF FLUCTUATION-DISSIPATION RELATIONS 

We consider a system of N discrete variables Ci (i.e. those entering models as Ising, Potts, Clock etc ...), generically 
called spins. Time t is discretized, namely i„ = nS, where n is an integer, and the time-step is 5 1/N. A 
configuration update is attempted at each time step. 

A. Transition probabilities 

Spin variables evolve in discrete time according to a generic Markov chain regulated by the transition probabilities 
w{a"\a' ,n) to go from a configuration a' to another a" in the n-th time-step. Transition probabilities obey the 
instantaneous detailed balance 

w(cr"|cr',n)exp[-/3-H(CT',n)] = w{a'\a" ,n) exp[~pnia" ,n)], (1) 

where 'H{a,n) is the (time dependent) Hamiltonian of the system. The diagonal terms w{a'\a\n) remain fixed by 
the normalization condition 

w{a'\a',n) = l-^w{a\a',n). (2) 

Restricting, for simplicity, to the case of single spin update, the form of the transition probabilities at time n is 

w{a"\a',n)^^J2'^,{<j"\a',n), (3) 

k 
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where Wk are the single-spin transition probabihties, namely a" and a' may differ only for the fc-th spin. 

The two-time conditional probability P{(y, n\a' , ni) to go from a' at time m to ct at time n can be expressed as 

F(a,n|a',m) = ^^ E (ala^-D, n) . . . u;„„ (a^^+D m) . (4) 

j„_i,...,j,„ <j("-i),...,o-(m + i) 

In the case of time independent w, the conditional probability is time translation invariant. For later use, we write 
this property as 

P((T,n|cr',TO+ 1) P(cr,n - l|cr',m). (5) 

Given two generic observables A{a) and B(a) (namely functions of a configuration of the system), from the knowl- 
edge of the conditional probability one can compute their correlation function 

= {A{n)B{m)) = ^ A{a)P{a, n\a' ,m)B{a')P{a' , m). (6) 

a, a' 



B. Relation between perturbed and unperturbed transition probabilities 



In the presence of an external perturbation hj{n) switched-on in the j-th site, the evolution is controlled by the 
Hamiltonian 'H(a,n) — JiQ^a) — ajhj{n). In the following we will always consider time- independent unperturbed 
transition probabilities and we will drop the time dependence in the unperturbed transition rates. The detailed 
balance condition ([1]) for the perturbed transition probabilities reads 



where, from now on, Wj and refer to unperturbed and perturbed transition probabilities, respectively. The most 
general form of u;^ obeying Eq. ^ is 



•j^{a"\a',n) = «;,((t'V', n)e^'''(")K— ^■)Af,(CT', a", n). 



(8) 



where Mj(a' ,a" ^n) is an ^.-dependent function symmetric with respect to the exchange of its arguments and such 
that Wj is a probability, namely positive and normalizable. To linear order in the external perturbation one has 



1 - ^hj{n){a'j - a'-) + mj{a" ,a')hj{n) 



1 — —hj{n){aj — a-j) + mj{a,a' ,n)hj(n) 



5.' 



(9) 



where 



Tnj{a" , a' , n) 



dMj{a",a',n) 



dhj (n) 



h=0 



(10) 



Let us comment on rrij: It is well known that the detailed balance condition leaves an arbitrariety on the choice of 
the transition probabilities, both for the unperturbed and the perturbed ones. Even for a fixed choice of unperturbed 
Wj, therefore, there is a family of different Wj obeying detailed balance, parametrized by rrij. 



C. Response function 



For a magnetic perturbing field hj{m) turned on the site j in the m-th time-step the impulsive response function 
Ri_j{n,m), describing the effect of the perturbation on the spin cTi at time n > m, is defined by 



J (n, m) 



1 d{<J,{n)) 
6 dhj (to) 



N 



h=0 



dhj (to) 



(11) 



h=0 
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where averages (. . .) are taken over thermal histories and the initial condition. 
From Eq.Q, one has 



^ „ dw'l{a"\(T') 

Rij {n, m 



P{a',m). 



(12) 



h=0 



The derivative of the with respect to the field can be easily obtained from Eq.Q 



dw'l{a"\a',n) 



dhj (to) 



h=0 



where 



and 



(13) 



9j (o-) Wj (cr|cr) = - ^ {a'\a) fj [a', a) 

g'^<7 



(14) 



(15) 



From Eqs. (|13ll4ll5p it is clear that the response function of Eq. ([T2|) cannot be straightforwardly interpreted as 
correlation functions. In order to do that, one would need the full transition probability /'(u, n\(j" ^ to + l)w{a"\'j' , m) 
connecting a' at time to to a" at time n, with w{a"\a' ,m) containing all the 'Wk{<y"\<y') according to Eq.Q, while 
in Eqs. (I13p only the one site 'Wj{(7"\(7') appears. A way out is to insert the missing 'w{a"\a' ,m) by writing 
NdWj- /dhjlh^Q = dw^ /dhj\}i^Q — u)(i(ln w'')/d/ij|^=0 5 as proposed in [l^, obtaining 



i?i J (n, m) 



E 

a, a' .a' 



a^P{a,n\a" ,m + l)w{a"\a') - ^ ' ' 



dhi 



P{a',m) 



(16) 



h=0 



However, let us notice that, although the response function is expressed in terms of the unperturbed dynamics, the 
function appearing on the r.h.s. of Eq. (jl6p is not in the form of a correlation function between observables according 
to the definition ([6]). This is because d{\nw'^) /dhj\fi=o depends on two configurations. 

Going back to Eq. (IT^ . in order to illustrate the CRT and the LCZ approaches, it is useful to write the response 
function as the sum of an off-diagonal contribution Di,j{n, m) and a diagonal contribution Dij{n, m) 



Rij{n,m) = Dij{n,m) + Di,j(n,m) 



(17) 



with 



and 



A:j(n,TO)-7V ^ (7,P{(7,n\a",m + l)wj{a"\a')fj{(7",a')[l-S,„^^,]P{a',m) (18) 

(7, a' ,(t" 



D,j{n,m) = N ^ a^P{a,n\a" ,m + l)wj {a"\cr') g-j {a') 6^" ,a'Picr' ,i 



(19) 



The above equations are exact and fully general. The next step is to express Di j and Di j in terms of correlation 
functions of observable quantities. This can be done in two different ways, leading to the CRT and LCZ results. We 
describe them separately below. 



D. CRT class algorithms 



Given the time interval (n,m), in numerical simulations one fixes a sequence I{n,m) of sites to be updated and 
then sums over different sequences. This corresponds to rewrite the conditional probability Q in the form 



P(ct, n\a' , to) — 



J\J'7i—rn Z ^ 



(20) 
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where the sum extends over ah iV" possible choices of X(n, m) in the interval [m, n]. Hence, {l/N)P{(j, n\a" , m + 
l)'Wj (ct"|(t') is the conditional probability restricted on the ensemble of trajectories satisfying the constraint /(m) — j, 
where /(m) is the particular site updated at time m in a given trajectory. This implies that Dij{n, m) can be written 
as the correlation {(yi{n)fj(m)5i(^rn).i) flip between Ui and fj, taking into account only trajectories where the j-th 
spin has been flipped at time m. Similarly, Di,j{n,m) is the correlation {ai{n)gj{'m)6j(rn),j)nofUp between ai and gj 
including only trajectories where flipping aj has been attempted at time m but rejected. Hence, the response function 
can be written as 



Rij{n,m) = N{ai{n)fj{m)Si(^^-)j)fUp + N{ai{n)gj{m)S^jn),j)nofiip- (21) 

This result is fully general. It holds irrespective of the nature of the discrete variables and of the form of the transition 
probabilities w and w'^. Notice that, because of the 6 function, on average only one out of N trajectories contributes 
to Rij- Therefore the overall factor N makes Rij well defined in the N oo limit. 

Chatelain Q and Ricci-Tersenghi [7] have considered the particular case of Ising spins interacting via the Hamilto- 
nian H{a, n) = — ai[Hj^ {a) + where Hj^ {a) = JJ2<i>j is the Weiss field (the sum runs over the spins 

interacting with ) , and of heat-bath transition probabilities 

K, n , exp[^(gf (a) + /^,-(m))a;.] 
"^-("'"'"^) = 2cosh[/3(gr(^) + /^.-(^))] - ^''^ 



This specific choice corresponds to 



/,(a,a')=/?(a,-(7f)=5,(^), (23) 



where aj ~ tanh(/3_ffj ), allowing to rewrite Rij{n,m.) in the more compact form 

i?,j(n, m) = N/3{a,in) [a,(m + 1) - (m)] <5z(™)j). (24) 

Here, since fj = gj, the distinction between {. . .) fup and {■■■)nofiip in Eq. (|2ip can be avoided. Although Ri^j in 
Eq. ([24l) is related to averages in the unperturbed dynamics, the ^/(m) j" acts like a projector on the restricted ensemble 
of phase-space trajectories including an attempted update of aj at time m. This is also the ensemble of trajectories 
that contributes to Rij{n,m) in standard numerical simulations where the perturbation is applied. The presence of 
the projector Sj(^jyij j makes necessary the knowledge of the sequences of updated spins restricting the applicability of 
this FDR to numerical simulations. This problem is bypassed in the LCZ algorithm, as shown below. 



E. LCZ algorithm 

In this section we re-derive the results of rcfs. fs', T^], originally obtained in a continuous time formalism, in the 
case of evolution in discrete time. Starting from the definition (I19p of Dij and using the time translation invariance 
property ([S]), one has J (n, to) = ^, criP(cr, n— 1|(t', m)_Bj((T')P((T', to) where i?j(iT) — {f3/2)wj{a\<T)gj{(T). Hence, 

A J {n, to) = (a, {n - 1)B^ (m)) . (25) 

We stress that, differently from the CRT scheme of Eq. ([M)) . the above form implies that no projection over a 
restricted ensemble of trajectories is present. 

We now turn to consider Dij. To begin with, taking advantage of the arbitrariness of TOj, let us consider the 
simplest choice rrij = in Eq. The effects of different choices of rrij will be considered in Sec. Ill Fl Then, from 
Eq. one has fj{(j", cr') = — (/3/2)(ct^ — a'-) and, since a' and a" may differ at most for the spin on site j, one can 
write 

^w, {a"\a') K - CT^) = ^ E (-'V') K (-'V') K - ' (26) 

k 

showing that Wj can be replaced with the full transition probability. Inserting into Eq. (|18p . Di,j{n, to) takes the form 

D^j{n,m) = ^{a^{n)Aaj{m)), (27) 
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where 

Aaj{m) = N[aj{m + 1) - crj(m)] (28) 

allows one to identify the discrete time derivative with respect to m of the autocorrelation function C{n,m) — 
{aj{n)aj{m)) in Eq. (P7)) . 
In conclusion, with the choice nij = made in Q, one has the relation 

= I [{<j,{n)Aa,{m)} - - l)B,(m))] , (29) 



with 



i3.(a) =^m,(aV)(a:-aO, (30) 



which is the form usually considered in the applications [1, [T^ US] • Notice that B depends on a single configuration 
and hence the term involving it in Eq. (j29|) is a correlation between observable quantities. 

As stressed previously, the above result, in addition to being general with respect to the form of the single spin flip 
unperturbed transition probabilities, holds true @ also for transition probabilities involving multiple-spin updates (as, 
for instance, Kawasaki spin-exchange). Extensions to the response of generic observables and to the case of transition 
probabilities that do not obey detailed balance are discussed in [ll|, [ij] and in [l^ , respectively. 

In Eq. (|29l) . at variance with the CRT result, no reference is made to the site I{m) to be updated at time to, 
and therefore there is no restriction on the ensemble of trajectories to be considered. The average over all possible 
choices of /(to) is, therefore, analytically performed. As it will be shown in Sec. Illll this makes the LCZ more efficient 
in numerical applications. More important, being an ordinary non equilibrium average, Eq. (j29p is well suited to 
standard analytical calculations and, in principle, to experiments. 

Finally, let us point out a property of correlations involving Bi that will be useful in the following. Given a generic 
observable 0{m) at a time m < rt — 1, from the definition (15111) one has 



{B,{n)0{m)) - (Aa,(n)0(TO)). (31) 



(B,(n)0(TO))= w,{d\a){d,~a{)P{a,nW,m)0{cj')P{o',m) (32) 



Indeed, 



and using Eq. (|26| one obtains Eq. (ISTI) . Eq. (|3T|) shows that in the mean Bi plays the role of the time derivative of 
a spin. 

F. Extra contributions related to rrij 7^ 

We now explore the consequences of a different choice of mj 7^ within the the LCZ scheme. Retaining the rrij 
contributions in Eq. Q, the response function can be written as 

i?jj(n, m) = Rf^f^in, to) -|- e,j (n, to) (33) 

with 

eij{n,m) — jS ^ (7iP{a,n\a" ,m + l)^Wj{a"\a')mj{a" ,a') [I — da",a'] 

a,(T" .cr' 

+ J2 m,{a",a')w, (a|a') <5,"..'}P(a', to). (34) 

Since the full transition probability cannot be reconstructed as in Eq. (I^Bl) . Sij can be identified as a correlation only 
in the restricted phase space of trajectories with the j-th spin updated at the time to as in the CRT scheme. The 
choice rrij = 0, therefore, has the advantage of avoiding this problem. 

It must be stressed that the formal manipulations leading to Eqs. (1^^ and p3l) are exact and hence they are 
identical if the same transition probabilities w^, or equivalently the same choice of M, are considered. In particular. 
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Eq. (p3)) contains the CRT relation of refs. [1, 0| as a particular case when heat-bath transition probability (|22l) . 
corresponding to ■mj{a",a') — — crj^(cr') is used. 

Therefore, let us compare the CRT and LCZ results in this case. Observing that, from Eq.Q, J^a^a' ('''1'^') = 
1 — Wj{(T'\a') and using Eq. ([5]), one has 

eij{n,m) = -iV(a-i(n)cr]^(m)(5/(„i) j) + {(Ti{n - l)(T]^(m)). (35) 

The first observation is that ei,j{n,m) = in equilibrium. Indeed, from the definition of crj^ one has 

{(Ti{n)<TY {'rn)5i{^m^ j) — {<Ti{n)aY {'rn + l)5/(m)j). Then, we use time reversal invariance to exchange the time ar- 
guments. The S function acting now at the larger time can be replaced by a factor 1/N representing the fraction 
of contributing trajectories. Finally, exchanging again the time arguments, one obtains N {ai{n)aY {'^)5i{m),j) = 
{ai{n — l)crJ^(TO)) and the right hand side of the above equation vanishes. Out of equilibrium this is no more true. 
However, it is generally expected that large-scale, long-time properties of scaling systems in the thermodynamic limit 
are not affected by the precise form of transition probabilities, provided detailed balance hold. The effect of different 
choices of m^, therefore, is expected to be negligible. Numerical simulations, presented in the next Section, confirm 
the expectation. 



III. VARIANCES 



The FDR's of CRT and LCZ have opened the way to numerical algorithms for the computation of the response 
function without applying the perturbation, the so called field-free methods. It was shown in Refs. 0, [HI that the 
calculation of the response function made via the CRT and LCZ algorithms is very precise and numerically efficient. 

In this Section we compute analytically the variances of the fluctuations of the response function obtained with 
the standard method (SM) where the perturbation is switched on, and with the two field-free methods. This task is 
carried out for Ising spins, the CRT method being valid only in this case. This allows us to compare the numerical 
efficiency of the different algorithms and to comment on the physical relevance of the variances, particularly in the 
context of systems with quenched disorder (see Section fill Bl) . 

Let us start by defining the fluctuating response function r^.j by 



and, therefore, its variance by 



where 



(n, m) = {vi^j (n, rn)) (36) 



Ag^ (n, to) = Rf} (n, to) - ^ (n, m) (37) 



RiAn,m) = {nj{n,m)nj{n,m)). (38) 



(2) 

We then focus on R^ J , computing it separately in the three methods. 
1. Standard method 

In the standard method one applies a sufficiently small magnetic field h at time to in the j-th site, and the 
response function is obtained by numerically implementing Eq. (jl2p where 



dw'}{a"\a') 



dhj 



w';ia"\a')-w,{a"\a') ^^^^ 



hi 

h=0 •' 



Since w^{a"\a') enters Eq. (fT2|) as the probability to flip aj at the time m, the above numerical derivative 
takes contribution different from zero only on the ensemble of trajectories were at that time the update of Uj 
is attempted. Then, imposing this restriction by means of the projector (5/(,„) and taking into account that 
((Ti(n)) = in the unperturbed dynamics, one obtains 

Rij[n,m) = N -^-^ — (40) 
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where the average is over the perturbed dynamics. We next observe that ((^/(m) j)^ = Si{m),j and that from Eq. 

(<5/(m)j) — 1/-^: since the 6 function cancels the sum over I{m) in P{a,n\a' ,n). From Eq. (HUl) one then 
obtains 

= (41) 

Notice that i?,- ^ diverges in the /i — ^ hmit. 

2. Ci?T relation 

Form Eq. (PT|) one obtains 

i?g(n,m) = iV^/?^ (/2 + l),a(™)) <5/(„.),,)^,,p + ^'Z?' (ffl ('T(m)) 5/(„o,.>„„/,,;p . (42) 

which holds true for any choice of rrij . Notice that the ensembles of trajectories contributing to the averages 
(. . ■) flip and (. . ■)nofiip SLie orthogonal. Therefore, no cross terms are present in Eq. (j42p . Restricting to the 
case of Ising spins and heat-bath transition probability, using Eq. (j23p one finds 

R^f'^Hn, m) = ^'N^[l - 2a,{m + l)af (m) + af (m)^] (43) 



3. LCZ relation 

In this case rij is given in Eq. (1^^ . and one has 

i?(f^^)(n,m) = (r,,,(n,m)r,,,(n,m)) = ^(Aa|(™)) + ^{B^im)) - ^(Aa,(™)i?,(m)). (44) 

Notice that (Acr|(m)) = 2A^^(1 — aj{m + l)crj(m)) = 2N'^{wj{(j' ^ cr|o')^7(m) j), where the last equality holds 
because only trajectories where at time m the j-th spin is flipped give a non-vanishing contribution. Hence, 
since the 5 function contributes on average only once every TV trajectories, (A(T|(m)) cx N . The second term on 
the r.h.s. of Eq. (j^J) docs not depend on N . Regarding the third term, reasoning along the same lines as for 
(AtT|(m)), from the definition ((28)) it follows that it is independent on N . Then, neglecting the last two terms 
in the large- limit one has 

i?(f^^)(n,m) = ^(Aa|(m)) = ^iV2(l - a,(m + l)a,(m)) (45) 



A. Comparison among variances 



As already mentioned, in the limit /i — >■ the standard method leads to a diverging variance. 

We now compare the variances of the two field-free methods using in both cases heat-bath unperturbed transition 
probabilities p2|). We first observe that (see Appendix I) 

1 - (a,(m + l)-7,(m)) - ^ [l - (a,(,7i)af (m))] (46) 

and therefore 

i?lf^^)(r.,m) = ^N{l-aY{m)a,{m)). (47) 
Next, for heat-bath transition probability, it can be shown (see Appendix I) that 

(a,(m + l)af (m)5,(„),,.) = ^(af(m)2). (48) 
Using the above result in Eq. P5|) . we get 

B^f'^^in, m) = P^N{1 ~ af {mf) , (49) 
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and finally, comparing with Ea. (l47l) . one obtains 

R^f''^\n,m) = 2i?f/^^)(n,m) + N^^ (frfM [a,(m) - (m)]) . (50) 

Recalling Eq. ([35]), and using Eq. (|48]), one can show that N(3'^ {a^ (m) [cri(m) - cr]^(m)]) = ej^j{m + l,m). As 
discussed in Sec. IIIF( this term is zero in equilibrium and one expects it to be negligible also out of equilibrium (this 
fact will be checked by numerical simulations in Sec. Iivp . Then one has 

-Rf/"^^^ (n, to) ~ 2i?f/^'') (n, to). (51) 

In order to compute the variances, according to Eq. ([57)1 the term ^ should be subtracted from R^i'^j ■ However, 

(2) 

these terms are negligible with respect to J in the thermodynamic limit being independent on N. Hence 

^iR,CB.T)^^^^^ - 2A(5'^^^)(n,TO). (52) 

The numerical evaluation of fluctuation of response functions confirms the above result, as it will discussed in 
Section HVl 

The origin of the factor 2 in the variances can be related to the different ways the term Dij of Eq. (fT9|) is treated 
in the CRT and LCZ methods and, in particular, to the presence of the (5- function in the CRT scheme (1^^ . Indeed, 
from Eq. (|2T|) one has that, in the CRT scheme, the fluctuating part of both Dij and Di j (corresponding to the two 
terms on the r.h.s) are non vanishing only once every N trajectories, and in this case their contribution is of order 
A^. Therefore, both the contributions to the variance coming from Di j and Dij are of order TV. Conversely, in the 
LCZ scheme, this is true only for the term Dij of Eq. (l27t . because Aaj in Eq. (|28|) is of order TV only once every 
TV trajectories when the j-th spin is flipped. Instead, from Eq. ([25]) one has that all trajectories provide a term of 
order one to Di j. So, the contribution to the variance associated with this term is of order one, and hence negligible. 
In conclusion, in the thermodynamic limit there are two terms contributing in the same way to the variance for the 
CRT algorithm whereas only one survives in that of LCZ. 



B. Integrated response function 

As already mentioned, the measurement of the impulsive response function R is numerically very demanding, so, 
in order to reduce the noise, usually the time integrated response function (dynamic susceptibility) is considered 

1 " 

Xi,j{n,m) = X! -Rij("'0 = (53) 

l—m 

where cc^j is the fluctuating part of Xi.j- In numerical simulations we focus on the equal site integrated re- 
sponse Xi,i- Taking advantage of space translation invariance, one usually computes the spatial average x{n^m) = 
(1/TV) X^iLi Xi^i{n,m) which fluctuates less than Xi.i. The variance of this quantity can be written in the form 

A^^'>in,m) = A[,^^(7i,m) + A(^)(n,TO), (54) 

where 

A[,'^)(n,TO) - - ^x'in,m) (55) 

1=1 

contains only equal sites terms, and 

Ai^\n,m) = -^"^{xt^iXj,]} ~ ^ ^ ^ x^{n,m) (56) 

is the contribution from different sites. 

In the case of simulations with the external field, the standard procedure consists in switching on a random 
perturbation during the interval [n, m]. One generally uses the bimodal distribution hihj = h?5i,j where the over-line 
indicates averages over the external perturbation. The integrated response function is then given by [l3| 

I N 

X{n, m) = ^ {(Ji)hhr (57) 

i=l 
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where Qh is the average in the presence of the perturbation. From the above equation and using the bimodal 
distribution of the external field one obtains 



and 



A^^^(n,m) = ^-lx'K™) (58) 
A('^)(n,m) = J;^^{siSj)hh^hj - ^ ^ ^ x^{n,m), (59) 

that can be also written as 

A(x)(n,TO)=^xg'^Km), (60) 

where 



- X^,i(n,m)xjj(n,m), (61) 

is a second order susceptibility. This quantity represents a tool for identifying cooperative effects in disordered systems, 
as it was proposed in [111, [13, and checked numerically in [IT, IS] . 

We then turn to consider the algorithms without the probing field. We first observe that the term is identical for 
all the algorithms and is always related to the non-linear susceptibility X^'/j Eq. (j60p . Indeed, from the definition 



(f56|) one has 

1 " " N - ^ 

Ai^Hn^m) = E E {r^AnJ)r,An,n) " ^^x'(",m). (62) 

i^j l—m I' —m 

The term {ri^i{n,l)rj,j{n,l)) in Eq. (|62|) can be written, using Eq. (IT2l) . as 



{ri,i{n,l)rjj{n,l)) = E criCrjP(cr, nja", Z + 1) 



dhi 

(7, a" ,(t' ,(T,a' 



P{a',l\aJ' + 1) 



P{d',l'). (63) 

h=0 



h=0 



The r.h.s. of this equation can be readily interpreted as the second order response Rf'-'^\n^ I, I') = ^g'}^.^(l"g'^^.(]!f 

leading to Ea.(l60|. On the other hand, the term Ao has different behaviors for the different algorithms. As already 
shown, Aq diverges as ft, — in the standard method. We now explicitly consider the term Aq in the CRT and LCZ 
algorithm. From the definition 

^ N n n ^ 

Ai''\n,m) = ^E E E {ndn,l)ndn,n) - ^x'(n,m). (64) 

i—l l—m I' —m 

As shown in Appendix II, {rij{n,l)ri_j{n,l')) = for / ^ I' and, therefore, only the terms with /' = I contribute in 
the double sum in this equation, yielding 

N n 

Ai^\n,m) = ^EE^^^'O- ]^X'(^.™)- (65) 

i—l I — 771 

In both the CRT and LCZ algorithms A^,^' is an increasing function of time roughly proportional to n — m, as already 
pointed out in 7]. This follows from substituting Eq. (|45l) into Eq. (|55|) and using Eq. p8)) . obtaining 

A('^)(n,m) = ^ E[l - (^^'^(0')] - -^X\n,m). (66) 

l—ra 

Since, at finite temperature, (o']^(l)^) is strictly less than one, the first term gives a contribution growing as n — m 
whereas x^ is always at most equal to /3^, and then sub-dominant at large times. Therefore, from the result (|52p at 
large times one has 

/\^^'^^^\n,m) = 2l\^^^^^^\n,m). (67) 
The numerical analysis supporting this result is presented in the following Section. 
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IV. NUMERICS 

In this Section we present the results of the numerical computation of the integrated response function Xi,i of Eq. 
(1551) using the SM with heat bath transition probabihties, and the CRT and LCZ methods. Details on the numerical 
implementation of the algorithms are given in Appendix III. 

On the basis of the analysis of Sec. |lTl since the LCZ algorithm corresponds to a different choice of Mj in Eq. ([8]) 
with respect to the other two, one may expect some differences in the results. However the numerical data presented 
below show that the response function computed with all the methods is the same within the numerical uncertainty. 
We then compute the variances of the response function, in order to check the analysis discussed in Sec. IIIII and to 
quantify the performance of the different methods. In the second part of this Section we also present data for the 
Fredrickson-Andersen model, showing that the LCZ field-free algorithm can be successfully applied also in this case. 

A. Ising and EA models 

We consider N — 100^ Ising spins interacting via the Hamiltonian 'H(cr) = ~ j) Jij^^i'^j, where the sum runs over 
nearest-neighbour spins on a three-dimensional cubic lattice with periodic boundary conditions, evolving according 
to heath-bath transition probabilities. The quantity B entering Eq. (|29)) takes the form 

B,{a)=a^-a,. (68) 

In particular we focus on the ferromagnetic Ising model (J,j — J ~ 1) and on the EA model ( = ±1 with equal 
probability). Temperature is measured in units of J. 

Let us start with the Ising model. In Fig. [1] the off equilibrium evolution of the system after a quench from infinite 
temperature to Tc = 4.5115 is considered. The susceptibility computed with the CRT and LCZ methods, and with the 
SM with h = 0.1 and h = 0.5 is plotted against n — m in the left panel. As it can be seen, the first three computations 
yield the same result with good accuracy. The SM with h = 0.5, instead, agrees with the other cases only up to 
n — m ~ 3 • lO''. This shows that non-linear effects become important from n — m ~ 3 • 10^ onwards. x(?t., m) grows 
monotonously to the equilibrium value Txeq = 1 as already found in previous studies [l^ . 

The right panel shows the variances A'^'<^^ Results with the SM, both with h = 0.1 and h — 0.5, give a time- 
independent value very weh consistent with A^^.^a/) ^ A^,^"^*^^ l/{Nh'^) (see Eq.[5Hl). This implies that Af'^"^*^^ 
is negligible and then Ag'^''^^^-' dominates the fiuctuation of x- With the CRT and LCZ algorithms one finds variances 
which are proportional to each other, within the numerical uncertainty, i.e. A^-^'*^^^' (n, to) ~ 2A^-^'^*-^'^^(n, to) in 
agreement with Eq. (j52p . and, as expected, they grow approximately linearly in n — to. According to the analysis of 
Sec. nil Bl this implies that Af'^-' is negligible with respect to Aq'^-'. Actually, this can be checked in Fig. [TJ where the 

term Aq'^'' alone is plotted, showing that it substantially coincides with the whole variance A(x) both for CRT and 
LCZ. 

Notice that, since A^^'^*^^ is constant while A^^'*-^^-^) and A^^'^*-^^^ grow in time, A*^^'^^^^ becomes smaller than 
the other two for large times, as it can be seen in the case h — 0.5 in Figure [TJ However, when this happens, the 
data obtained with the SM are already affected by large non-linear effects, as it is evident from the difference between 
the curves corresponding to ft. = 0.1 and h — 0.5, in the left panel. Hence, the signal to noise ratio in the field-free 
methods is better than the one in the SM, provided that one works in the linear regime, as it is shown in the right 
panel. Moreover, the factor 2 between the variances implies that, in order to have a certain signal to noise ratio, 
simulations performed with the CRT algorithm require a larger statistics (by a factor \/2) than those based on the 
LCZ method. 

Fig. [5] displays the behavior of the same quantities as in Fig. [T]in the case of a quench to T = 3, below T^. The 
behavior of the susceptibility is now characterized by a maximum around n — to ~ to, due to the interplay between the 
response of single interfaces and the reduction of their number [l7j , as domain coarsening goes on. For what concerns 
the comparison between the variances, we first observe that in the SM the variance is almost time-independent around 
a value in good agreement with l/{Nh'^). With the field-free methods, AQ'^-'(n, to) still grows linearly with n — m 
and the relation (p7)) is very well verified. Nevertheless, at variance with the quench to Tc, Aq'^-' does not represent 
the whole variance A^-^-* and the contribution Ar'^'' is not negligible. As already mentioned, we expect Ar"*^-* to be the 
same for all the methods. This is shown in Fig. [Sj where Ar"^"* = A'^'<^) — Aq'^'' is plotted against n — m. The behavior 
of A^'^'' is consistent with the power law Ai^\n,m) oc (n — to)° '''. Therefore, since the growth of Ar"*^-* is slower than 
that observed for Ag'''', the proportionality A^^^'^^-'") ~ 2A(^'-^'^^-' between the whole variances is expected to hold 
at times larger than those accessed in the simulation. 
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n-m n-m 



FIG. 1: (Color online). Integrated auto-response function x(Wi ^) (left panel) and its variance (right panel) in the 3-dimensional 
Ising model quenched from T = oo to Tc = 4.5115. Left panel: Different curves correspond to computations performed with 
the SM (using two different values of h) and with the field- free algorithms of CRT and LCZ, as indicated in the key. In the 
right panel the behavior of the variances of the response function computed with the different methods is shown (continuous 
lines with heavy symbols, see key). The terms Aq'^''^^"^' and Aq'^'^'^'^' are also plotted (dashed lines). 



The observed difference in the behavior of A'*: in the quench to and to below Tc can be attributed to different 
structure of domains. In the quench to below Tc, domains are compact and inside domains {a^ (m)'^) ~ m^^ where 
nieq is the equihbrium value of the magnetization in the ordered phase. Therefore, indicating with p{t) the defect 
density and taking into account that at interfaces (cr]^(m)^) ~ one roughly expects {a'^ {m)'^) ~ fnlqi^ — p{t))- 
Substituting this result into Eq. (|66p one obtains asymptotically, when p{t) <C 1, 

A[,'^^(n,m) ~ (l-m^J (n-m). (69) 

Since in the ordered phase nieq is very close to one, this implies that Aq"^' linearly grows in time, as indeed observed, 
but with a very small prefactor (the smaller the lower is the temperature) . This makes Aq'^^ comparable or even sub- 
dominant at small time differences n — m with respect to A^"^^. Nevertheless, since A^'^'' grows with a smaller exponent 
(~ 0.5) fluctuations are always dominated by Aq'^-' at large times. Conversely, in the case of quenches to T = Tc, 

rrieq = and Ag'^'' is the dominant contribution even for small n — m. 

In Fig.|4]we show the integrated response function and its variance in the EA model quenched to (left) and below 
Tc (right). Here we take for Tc the value obtained in ref. ([1^). AH the methods yield the same result except the SM 
with the largest value of h which, for large values ofn — m {n — m^ 6 - 10® for T — Tc (left) or n — m ^ 10^ for T = 1 
(right)) is affected by non-linear effects. The variances show the behavior similar to the one already discussed for the 
ferromagnetic Ising model. A notable difference is that, not only in the critical quench, but also in the sub-critical case, 
one has A^'^^ ~ ^o'^^ indicating that aI'^'' is always a sub-dominant contribution. We have explicitly checked that 
this happens also in quenches to lower temperatures (T = 0.5, T = 0.2). This possibly indicates that the argument 
developed above for ferromagnets cannot be straightforwardly extended to the low temperature phase of spin-glasses 
and that second order susceptibilities (or variances) may be used as efficient tools to characterize this difference. 

B. Fredrickson- Andersen model 

As stated in [8] and discussed further in [ll|, [l^ and in Sec. |ll] of this Article, the FDR of LCZ can be derived 
in full generality in the context of Markovian dynamics, regardless of the model Hamiltonian and of the choice of 
transition probability. However, this claim has been questioned in 14] where the applicability of the LCZ algorithm 
to the FA model was doubted. However, the derivation of the LCZ relation as given above, or with different analytical 
techniques in Refs. [ill [13) [HI, shows clearly its general character where no particular assumptions on the specific 
system nor on the form of the perturbation are made. In order to illustrate this issue also numerically, we have carried 
out the computation the integrated response function with the LCZ algorithm in the one-dimensional FA model. 
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FIG. 2; (Color online). As in Figure [T] for the 3-dimensional Ising model quenched from infinite temperature to T = 3 < Tc. 
Left Panel: Different algorithms yield the same result except the SM with the largest value of h which, is always affected by 
non-linear effects. 




FIG. 3: (Color online). Difference A''^' (n, m) — Aq (n, m) for the three algorithms in the Ising model quenched to T = 3. The 
dashed orange line has slope 0.5. 



The Hamiltonian of the FA model reads H{(t) ~ J^i "'ij where the ai are bimodal variables taking the values (0, 1), 
with (Ti = 1 for a mobile fluid region and di = an for an immobile one. Spins evolve according to transition 
probabilities obeying detailed balance, whose off-diagonal elements are 

w,{a'\a) = [e(l - a,) + (1 - e)a,]A,(a), (70) 

where e — 1/ (l-t-e^^^) is the equilibrium density and \i{a) — Ui-i +ai+i is a kinetic constraint that preserves detailed 
balance, due to the independence from (7^. The integrated response function has been computed using both the SM 
and the LCZ algorithm. In the SM, the effect of the external field amounts to replace e with e'' = 1/(1 -I- e'^"''')/-^), 
whereas the quantity Bi defined in Eq. ([50)) and entering the LCZ relation is given by 

B,{<j) = [e(l - a,) + (1 - e)<j,]X,{a)[l - 2<j,]. (71) 

In Fig. [5] the comparison between the data obtained with the two methods and for three different temperatures, is 
presented. The value of h used with the SM was checked to be in the linear regime. In particular we found that keeping 
the ratio h/T = 0.01 constant satisfies the linearity requirement. The agreement between the two computations is 
excellent in the whole time range and for all the temperatures considered. Data for T = 1 converge to the equilibrium 
value whereas results for T = 0.5 and T — 0.2 exhibit a non-monotonic behavior already observed in one-dimensional 
kinetic constrained models [22*1. At low temperatures the equilibrium value is reached only asymptotically. 




FIG. 4: (Color online). The integrated auto-response function x{n,m) (insets) and its variance (main) in the 3-dimensional EA 
model quenched from T = oo to Tc — 1-2 (left panel) or to T = 1 < Tc. Insets: Different curves correspond to computations 
performed with the SM and two different values of h and with the field-free algorithms of CRT and LCZ, as indicated in the 
key. In the main part of the two panels the behavior of the variances of the response computed with the different methods is 
shown (continuous lines with heavy symbols, see key). The terms Aq'^''^^"^' and Ag'^'^'^'^' are also plotted (dashed lines). 




FIG. 5: (Color online). Results of the numerical computation of the integrated auto-response function xi^jm) for the one- 
dimensional FA model with m = 10'''. The response functions computed with the SM (red filled symbols) and with the LCZ 
algorithm (continuous black lines), agree for the three temperatures. In the SM the external field is h = O.OIT, satisfying 
linearity. The equilibrium value (green empty symbols) is also plotted for each temperature. 



V. CONCLUSIONS AND PERSPECTIVES 



In this paper we have compared the FDR derived in a series of papers by CRT [1,0] and LCZ [8]. First, by re-derivinff 
them in a unified formalism we have pointed out that the distinction between these FDR is due to a different choice 
of the perturbed transition probabilities. Actually, even restricting to systems where detailed balance is obeyed, for a 
given choice of unperturbed transition probabilities an arbitrarily remains on the form of the perturbed ones, which is 
parametrized by the function rrij defined in Eq. PT7)) . In the case of Ising spins, selecting the value of rrij corresponding 
to heat bath transition probabilities leads to the CRT relation (l24l) . This FDR relates Ri,j to a correlation function 
involving a (5-function which weights only a subset of the whole ensemble of unperturbed trajectories. For this reason, 
the FDR of CRT, besides being limited to Ising spins with heat bath transition probabilities, can only be used in 
numerical simulations. On the other hand, making the choice mj = of LCZ allows a further mathematical treatment 
leading to the FDR p9| where the response function is related to standard unperturbed correlations. This makes 



15 



this relation basically different from those obtained in fld\ (Eq. p6|) and in other approaches (i.e. in [^Q), where 
the response function cannot be expressed in terms of correlations of observables. This makes the applicability of the 
LCZ relation in principle not restricted to simulations. Furthermore, this FDR has a larger degree of applicability 
with respect to the CRT, since it is not restricted to Ising spins nor to heat bath unperturbed transition probabilities. 

In the second part of this Article, we have studied the efficiency of the CRT and LCZ field-free methods. In order 
to do that, we evaluate analytically the variances of the response function obtained with the SM or with the field- free 
methods. It turns out that, as far as the autoresponse function is considered, field-free methods are by far more 
efficient than the SM. This combines with the advantage of having linearity (h — 0) built in. Moreover, we found 
that the LCZ algorithm is slightly more efficient (by a factor than the method of CRT. 

We conclude by pointing out that the results contained in this paper are not restricted to the framework of the 
efficient computation of the response function. Indeed, we mention that the study of the variances is closely connected 
to the issue of characterizing the fluctuation of two-time quantities in aging systems, a problem which has received a 
good deal of attention recently Q. Moreover, as discussed at the end of Sec. Illli the variance of the response function 
is also related to the second order susceptibility introduced in [ll|, [ill study of cooperativity. 
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VI. APPENDIX I 

We first prove Ea. (|46)) 

1 - {a,{m + l)a,{m)) = ^ [l - {<J^ {m)a,{m))] . (72) 
To do that, we first observe that from Eq.(l3]) 

(l-a,(m + lK(m))^l ^ [1 - afa^, (ci'V) ^^(^', (73) 

a" ,(t' 

Only the term Wi {ct"\(j') in the transition probability between the time m and m -\- 1 contributes. Hence the two 
configurations cr" and cr' differ only for the spin on the i — th site and sum on a" reduces to the sum on the two 
possible values icr,'. Then, using the Heath bath form for the transition probabilities, 

w,{a"W)^\{l + <j';aY{a')) (74) 

one obtains Eq. ([7^ . 

Next we prove Eq. . Let us first proof that 

([a,(m + 1) - (m)]af (m)^,(„),,) = 0. (75) 

From the definition one has 

([a,(m + 1) - (m)] {m)5i(^^)^,) = ^ - (a')] w, (a'V) 5/(™),,P(a', m). (76) 

cr" ,a' 

Because of the delta function, one has that the configurations a' and a" can differ only for the spin in the j-th site. 
This implies that gj^(cr') — {o'") and that the sum on a" reduces to the sum on the two possible values ±ct^. The 
Heath-Bath form (|74|) for the w then gives 

[a;'-af(a')]u;,(a'V')=0, (77) 
which implies Eq. ([75)) . We next observe that from Ea.d^O]) 

I{m) a' ,a 

Combining the above result with Eq. ([751) one recovers Ea. (H51) . 
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VII. APPENDIX II 



Here we prove that R\J (n, m, m') — {rij(n, m)rij(n, m')) is identically zero for all the algorithras if m 7^ m'. 
1. Standard method 



In the SM one applies a random magnetic field /i^, with hi = and hihj = Sij, and the response is obtained as 



Ri,j[n,m)^N ^-^ . (79) 



Therefore 

/l2 



(2) N 

R\ ' {n, m,m') = -rnSm,m' (80) 



2. Ci?T relation 

In this case r^j- can be directly read off from Eq. (|2ip . From Eq. (|T5|) one has 

^ u;,(a'V')/.K, - <5.",.']P(fT', to|(t, to') = - 51 M^'W')9^{^')Sa".a'PicT', m\a, to') (81) 

a" ,<7' a' 

that gives 

(/,(m)0(m')) + {g^{m)0{m'))„ofl^p = 0, (82) 
for every generic observable O computed at the shorter time to'. Therefore, 

i?g(n,TO,TO') =0 (83) 

for m' < TO — 1 . 

3. LCZ relation 

Reading r^j- from Eq. (f29| . using the property (|3T|) one obtains that {[/S.ai{m) — i?i(TO)]0(m')) = 0, for every 
generic observable O computed at the shorter time m' . Therefore one arrives again at Eq. (|83p for to' < to — 1 . 



VIII. APPENDIX III: NUMERICAL IMPEMENTATION OF ALGHORITHMS FOR THE 
COMPUTATION OF INTEGRATED RESPONSE FUNCTIONS 



1. Standard method 

For each realization of the dynamical trajectory, a random magnetic field hi is assigned to each site, hi is 
usually chosen from a bimodal distribution hi = ±h. The evolution is then controlled by unperturbed transition 
probabilities until the time to and by the perturbed ones given in Eq.® for later times. At the time n, the 
integrated response function is computed according to Eg. dSTl) . 

2. CRT relation 

The integrated response function can be obtained from the space and time integral of Eq. (|24|) 

X^^^^Hn,m) = |,^(a.(n)A,(n,TO),) (84) 

i 

with 

A.(n,TO) = ^ h(/ + l)-ar(0]. (85) 

l—m 

In each realization of the dynamics, the quantity A^ is initially set to zero on each site. For all the timesteps 
I > TO, Aj is updated, via the relation 

A, = A,+a,(Z + l)-c7f (0, (86) 

only on the site j = I{1) where the flip of the spin has been attempted. Elsewhere A^ is left unchanged. At each 
time n>m, x(''T-,™) is then computed according to Eq. (j84p . 
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3. LCZ relation 



From Eq. (P^l one immediately obtains 



X 



2A^ 



[{(T.,{n)ai(n)) - {ai{n)ai{m)) + (CTi(n- l)Ai{n,m))] 



(87) 



with 



^,(n,TO) = -^B,(/). 



(88) 



l—m 



The quantity {fji{n)ai{m)) is the usual two time correlations function. Concerning the evaluation of Ai^ the 
basic observation is that, according to Eg. (1681) . Bi only depends on the spin ai and on the spins interacting with 
it. In particular, for the models considered in this paper, Bi depends on and on its nearest-neighboring spins. 
The evaluation of Ai proceeds as follows. At the time m, Bi is evaluated on each site according to Eg. (1551) . Ai 
is set to zero and U is set to m. li represents the time where the last evaluation of Bi has been performed. Bi 
is then left unchanged on all sites until a spin flip occurs. If the spin flip occurs on site I{1) at time I, Aj is 
updated as 



Ij is then set to / and the new value of Bj is evaluated. The same procedure is repeated for all the spins 
interacting with cr/(;), leaving unchanged Ai, Bi elsewhere. In the end Ai{n, m) is obtained via the relation 



and the integrated response function is computed according to Eq. (j87p . 
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